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Abstract 



We study lead nanodroplets containing 147 to 1415 atoms at temperatures 
ranging from the bulk melting point up to the beginning of the evaporation 
regime. Molecular dynamics simulation and an embedded atom potential are 
used. The structures, total energies, and mobility of atoms in the clusters are 
analyzed. We found that the liquid cluster of 147 atoms shows a pronounced 
tendency to form non-spherical shapes, and sometimes separates into two 
droplets. Bigger clusters disassemble by evaporation of monomers. We also 
explore shape oscillations of these nanodroplets using the nuclear liquid drop 
model. 
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I. INTRODUCTION 



Clusters assembled from metal particles are scientifically interesting and technologically 
important JT|. Physical phenomena in clusters of atoms and molecules do not necessarily 
have the same features as in bulk phases. Melting of metal clusters was already studied by 
computer simulations and in experiments pTOf] . In contrast, the studies of clusters 

at even higher temperatures, i.e., approaching the liquid- vapour transition are scarce. Ex- 
ceptions are early Monte Carlo and Molecular Dynamics (MD) simulations for very small 



Lennard- Jones clusters that extended to the temperature region of the vapour phase fll]--|14 



Related simulations were also performed in the field of nuclear fragmentation |]T3j|, and re- 



cently in the study of fragmentation for argon clusters fll6 |. Understanding the behavior of 



nanoparticles at such high temperatures may also shed light onto the process of achieving 



plasma and controlled fusion in a gas of clusters |L7[] . 

Melting of lead clusters of sizes ranging from a few to about 100 nm were investigated 
by high-sensitivity optical reflectance and dark-field electron microscopy JTU| . Spherical and 
non- spherical nanometric Pb inclusions in inert amorphous matrices of SiO and AI2O3 were 
studied and it was found that melting is enhanced by surface and curvature effects. These 
studies have been carried out up to the bulk melting point of lead: T s i = 600.56 K, while 
the liquid- vapour transition temperature for bulk lead is: Ti v = 2033.16 K. 

In this letter we present a MD simulation study for liquid magic number Pb clusters 
containing 147, 309, 561, and 1451 atoms. The high-temperature structures, total energies, 
the mean-square displacements, and diffusion of atoms are analyzed. A calculation of shape 
oscillations for lead nanodroplets in the liquid drop model has been carried out. 



II. METHOD AND COMPUTATIONAL DETAILS 

Properties of metals in bulk, surfaces and clusters are unlikely to be described well within 
the MD method without a many-body potential |fT8 |. We have used such a classical many- 



body potential for lead ||T9"|-f21|j . The set of parameters in the potential was determined by 



fitting to several measured properties of metal |22[| . This potential was already used in MD 
simulations of the structural and vibrational properties of lead clusters [|19|,|2~0| . The bulk 



melting point for the potential is T s i = 618 ± 4 K [pi]] . This value was precisely determined 
by simulating coexisting liquid and solid phases under constant energy. For the present 
investigation of the high-temperature behaviour of nanoparticles we have estimated the 
bulk liquid-vapour transition temperature as T[ v ~ 2050 K. The precise calculation of the 
liquid- vapour transition temperature for a given potential is in principle and computationally 
more involved (in comparison with the calculation of the bulk melting point) and will be 
discussed elsewhere. 

In this simulation lead atoms were at T = K distributed as regular Mackay icosaheadra 
of 147, 309, 561, and 1415 particles. In reference jn|, the authors find using the same 



potential that at T = face-centered cubic morphologies are favoured over icosahedra. This 
is the result of the low surface energy anisotropy between Pb(lll) and Pb(100), and of 
the high tensile surface stress which inhibits icosahedral structure. However, at T > the 
stress decreases and icosahedra become more favoured. In experiments with argon clusters 
and in simulations for clusters described by pair potentials (such as Lennard- Jones and 
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Morse), there is now clear evidence that for less than 1500 atoms icosahedra are more stable 
than fee structures [ft4| ) P5|l - The situation is less clear for metallic clusters. The majority 
of experimental and theoretical studies find that icosahedral clusters are also energetically 
preferred for metallic clusters of less than 1500 atoms, even at T = [26j,4|. Actually, the 
morphology of the energetically preferred structure for solid clusters at low temperature 
is not important for our study of liquid nanoparticles at high temperatures. Due to our 
wish to compare the results of simulations with future experimental studies of melting and 
evaporation we have chosen icosahedra as starting point in our simulation. As in other MD 
simulations of metal clusters [@-[7J] no boundary conditions have been employed. A time step 
of 7.33 x 10~ 14 s was used. MD simulations with damped dynamics were first performed 
to extract the kinetic energy until the clusters reached the structure corresponding to a 
local energy minimum. Then the clusters were heated up to temperatures (T < 2000 K) 
where evaporation of atoms starts within a typical time scale of classical MD simulations. 
The embedded atom potentials are known to show realistic evaporation, i.e., to be much 
less prone to enhanced evaporation during MD simulation in comparison with the pair 
potentials. First the "heating" runs at constant temperature were carried out for 5 x 10 4 
time steps. Because of the difficulties in achieving equilibrium in the liquid temperature 
region, rather long runs of 10 6 time steps (i.e., 7.3 ns) were used to calculate the cluster 
properties. The temperature was controlled by rescaling the particle velocities at each time 
step. Note that the potential used here (and other embedded atom potentials) are not 
applicable at extremely high temperatures, for example at T > 2000 K and approaching 
plasma conditions. These potentials are not fitted for this temperature region and the 
resulting low atomic coordination. 



III. RESULTS AND DISCUSSION 

A. Structures and mobility 

The equilibrium shape of a macroscopic liquid metal drop is a sphere because of the 
negligible anisotropy of the liquid-vapour interface free energy in comparison with the one 
of the solid-vapour interface. We have found that on the simulation time scale the liquid 
clusters exhibit substantial shape fluctuations at high temperatures. This tendency towards 
transient formation of non-spherical shapes decreases with increasing cluster size, but it is 
still present even for the largest cluster of 1415 atoms studied here. We may mention that 
non-spherical shapes were also found in experimental studies of nanometric lead inclusions 
T0|] • A snapshot of atomic positions producing a very elongated lead droplet is presented 



in Fig. |TJ (a). In the further time evolution this cluster separates into two smaller droplets 
as shown in Fig. [I] (b). The colour of the atoms reveals that the cluster emits a small 
subcluster. The process presented in Fig. |TJ is the only occurrence we found (within our 
long, but limited simulation time) for a droplet separation into two parts, although several 
other extremely elongated shapes similar to Fig. [IJ (a) were observed, especially for the 
smallest cluster and above T > 1000 K. For most liquid droplets a typical process at these 
temperatures is single atom evaporation. A similar tendency of less intense evaporation for 
small as compared to large nanoparticles was recently found in a MD study of melting of 
argon clusters f27|. Another rare event - a dimer evaporation was also detected. We have 
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found that the number of evaporated atoms strongly depends on the length of time evolution 
during the MD simulation. For example, a liquid nanoparticle of 1415 atoms at T = 1600 
K completely vaporizes after 3 x 10 5 time steps, whereas after 10 5 time steps only 4 atoms 
evaporate. For an almost immediate (within ~ 10 ps) complete "explosive" disintegration 
of lead nanoparticles of these sizes much higher temperatures of the order of 5000 K are 
needed. However, as pointed out above, the reliability of the potential used is questionable 
in this temperature region. 

In most bulk liquids the coordination number is ~ 8 — 11, i.e., it is smaller than the 
coordination number of, for example, 12 which is the property of close-packed bulk solid fee 
metals |[28| . In addition, as a result of surface effects, the average coordination number in 
clusters is always smaller than in the bulk. The average coordination numbers for all liquid 
nanoparticles studied here at several temperatures are given in Table |. Our calculation of 
the coordination numbers (and other cluster properties which we present below) extends 
to the beginning of the evaporation regime, i.e., up to T = 1200 K. We found that at 
T = 1000 K and T = 1200 K, after several nanoseconds, one or two atoms evaporated from 
the liquid clusters of 309, 561, and 1415 atoms. After such an evaporation process, the 
calculations were continued only for the remaining liquid cluster. From Table p] one can see 
how the average coordination number roughly increases with cluster size and decreases with 
temperature. 

Steps in the caloric curve are a signature of melting [|l]-[7j . We also studied the potential 
energy per atom, E, as a function of temperature from melting up to the evaporation 
regime. It was found that the E = f(T) function is linear for all studied nanoparticles in 
that temperature region. 

Atomic configurations presented in Fig. [I] show the mobility of the atoms on a qualitative 
level. The mean-square displacements for the clusters after several nanoseconds of simulation 
and at several temperatures are shown in Fig. At lower temperatures ( T < 1000 K), the 
mean-square displacements rise approximately linearly with time. Such a behaviour is typical 
for liquid diffusion. Sharp increases at higher temperatures are the results of evaporation, 
or the onset of evaporation. As pointed out above, the atom which has evaporated (and 
moved from the surface of the nanodroplet at a distance larger than the cutoff radius of 
the potential, i.e., > 5.6 A) is excluded from further consideration. Note that as a result 
of evaporation from the cluster of 561 atoms the mean-square displacements are bigger at 
1000 K (560 atoms in the liquid cluster), than at 1200 K (559 atoms). 

The diffusion coefficients of bulk solid and bulk liquid phases are orders of magnitude 
different (10 -9 cm 2 s _1 vs 10 -5 cm 2 s~ 1 ). The temperature dependence of the diffusion coef- 
ficient of lead clusters is shown in Fig. The values are typical for diffusion in the liquid 
phase. As for the mean-square displacements, a jump of the diffusion coefficient is the result 
of the onset of evaporation. 

B. Shape oscillations of liquid lead clusters 

It is interesting to estimate the influence of the size of a lead nanodroplet on its shape 
oscillations. Many phenomena in cluster physics can be analyzed using the tools devel- 
oped in nuclear physics p9[ . For nuclei some phenomena, for example magic numbers, can 
only be studied within independent particle models. Other phenomena, for example binding 
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energies, vibrational properties, and fission of heavy nuclei, are studied using collective (con- 
tinuum) models, such as the liquid drop model. Experimental and theoretical investigations 
have shown that the same two types of models can be applied for metallic clusters [E9| . 



The fission measurements for small gold clusters were interpreted using the nuclear liquid 
drop model and the results indicated that the vibrational properties of metallic clusters, 



nuclei, and liquid drops are similar [p0|| . A common approach to the analysis of vibrational 
properties of clusters is either by direct diagonalization of the dynamical matrix, or by cal- 
culating the Fourier transform of auto-correlation functions in a MD simulation. For liquid 
clusters bigger than 100 atoms as studied in this work, these traditional approaches give a 



broad, structureless line in the frequency spectrum pl| , p0| . Instead, we analyze the shape 
oscillations of liquid lead nanoparticles using the nuclear liquid drop model H32]] . In this 
model surface and volume modes are analyzed separately and the results are expressed by 
macroscopic properties of the metal. Three kinds of shape oscillations are possible: surface 
oscillations, modes of vibration involving compression, and coupled oscillations. 

The distance of the surface of a slightly deformed liquid drop from the origin, R(9,(f>), 
can be expanded as 

m <P) = Ro ^ + E axjGpie, <t>) j , (l) 

where i? is the equilibrium radius, are normal coordinates, and 1a m (6', <fr) are the spher- 
ical harmonics ||32|| . The equations of oscillatory motion are 

&xn + Lo(s) 2 x a XlM = 0. (2) 
Eigenfrequencies of the surface oscillations are given by 

= o\(\- 1) (A + 2) ^ (3) 

PO-Ko 

where a is the surface tension and po is the density In Eqs. (l)-(3), the value A = 
represents a compression without change of shape, whereas A = 1 describes a translation of 
a droplet as whole. This can be seen by considering Eq. (1) and respectively, the volume of 
an incompressible drop and a small translation of the center of mass along the axis 9 = 0. 
Therefore, the lowest order surface oscillations are quadrupole modes A = 2. 

Compressional oscillations in liquid droplets also exist. Small deviations of the density 
p(f) which are solutions of the linearized hydrodynamical equations that lead to the wave 
equation are given by 

Sp(f, t) = p ix{k nX r)Y^(e, (j>)a n ^(t), (4) 

where j\ is a spherical Bessel function. If the boundary condition is such that Sp(f,t) 
vanishes at r = Rq, then the eigenvalue equation is 

Jx(k nX R ) = 0. (5) 

The compressional eigenfrequencies are given by 
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w(c) n A = U c k nX , 



(6) 



where u c is the sound velocity in the liquid. 

The ratio of the eigenfrequencies of surface and compressional modes (Eqs. (3) and (6)) 
slowly decreases as the number of atoms in the cluster increases. In general, if the liquid 
cluster (or nucleus) is small, it is important to consider coupling of the surface and com- 
pressional modes. These coupled modes are solutions of the wave equation with a modified 
boundary condition such that the pressure induced by the shape oscillations is in equilib- 
rium with the pressure induced by the compressional oscillations p2| . The corresponding 
eigenvalue equation of the coupled modes is 



1 



j\(k nX R 



o 



d_ 

Or 



3\{K\r) 



=Ra 



Ro u(s) 



(7) 



A 



The surface, compressional, and coupled frequencies were calculated using the experi 
mental data for bulk lead: u s = 1800 ms^ 1 [Q, p " 1 nr;_v '•"•"~ ! 

Nm- 1 



10678 kgm~ 3 @], and a = 0.453 
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The average radii of approximately spherical clusters were taken from config- 
urations obtained in MD simulations at 600 K. The calculation shows that for these lead 
nanodroplets the frequencies of coupled modes are approximately equal to the frequencies of 
surface modes. The size dependence of the surface and compressional frequencies is presented 
in Fig. 4. Frequencies of compressional modes are of course bigger than the corresponding 
frequencies of surface modes. The frequencies increase with increasing value of A and with 
decreasing size of the nanodroplets. This agrees with our results of MD simulations where 
the most pronounced shape oscillations were found for the smallest liquid cluster of 147 
atoms. 

These results provide for the first time an insight into the behavior of neutral liquid 
metal clusters at high temperatures, up to the beginning of the evaporation regime. Exper- 
imental investigations of clusters in this temperature range, either using optical, diffraction 
techniques, and mass spectroscopy []]],§- nj, or cluster fragmentation properties |T6| , p!7|j35 
are desirable. 
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TABLES 



TABLE I. The average coordination numbers for liquid lead clusters of a given total number 
iV of particles, calculated after 7.3 ns of simulation. If evaporated atoms exist their number is 
given in parentheses; they were excluded in the calculation of the average coordination number 
and other cluster properties. 



Temperature (K) N = 147 N = 309 N = 561 N = 1415 

600 8^99 9W 9^34 9/72 

800 8.75 8.84 8.93 9.34 

1000 8.44 8.43 8.94(1) 9.04 

1200 8.41 8.13(1) 8.64(2) 8.88(2) 
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FIGURES 



FIG. 1. Atomic positions showing: (a) non-spherical shape of the liquid cluster of 147 atoms 
at T = 1300 K after 6.96 ns of time evolution, (b) the same sample after 7.32 ns. The relative 
darkness of an atom is proportional to its square displacement during the run. Black particles are 
those with the maximum of mobility. 

FIG. 2. Mean-square displacements for liquid lead nanoparticles at several temperatures: (a) 
147, (b) 309, (c) 561, (d) 1415 particles. 

FIG. 3. Diffusion coefficients for liquid lead clusters as a function of the temperature. 

FIG. 4. Eigenfrequencies of lead nanodroplets in the liquid drop model as a function size: (a) 
surface oscillations, (b) compressional oscillations (n = 1). 
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